#include "PoissonSolver.h"
#include "Grid.h"
#include "BoundaryConditions.h"
#include <iostream>
#include <fstream>
#include <nlohmann/json.hpp>  // JSON library
#include "exprtk.hpp"  // ExprTk 头文件
#include <cmath>
#include <vector>

using json = nlohmann::json;

// 解析字符串表达式并返回函数对象
std::function<double(double, double)> parseExpression(const std::string& exprStr) {
    return [exprStr](double x, double y) {
        // 创建符号表并添加变量
        exprtk::symbol_table<double> symbolTable;
        symbolTable.add_variable("x", x);
        symbolTable.add_variable("y", y);

        // 创建表达式对象
        exprtk::expression<double> expression;
        expression.register_symbol_table(symbolTable);

        // 解析表达式
        exprtk::parser<double> parser;
        if (!parser.compile(exprStr, expression)) {
            throw std::runtime_error("Failed to parse expression: " + exprStr);
        }

        // 返回计算结果
        return expression.value();
    };
}

// 计算误差的 1-范数、2-范数和∞-范数
void computeNorms(const Eigen::VectorXd& numericalSolution, const std::function<double(double, double)>& realFunction,
                  const Grid& grid, double& L1, double& L2, double& Linf) {
    int N = grid.getN();
    double sumL1 = 0.0;
    double sumL2 = 0.0;
    Linf = 0.0;

    for (int i = 0; i < N; ++i) {
        double x = grid.getX(i);
        double y = grid.getY(i);

        // 真实解
        double u_exact = realFunction(x, y);
        // 数值解
        double u_numerical = numericalSolution[i];

        // 误差
        double error = std::abs(u_numerical - u_exact);

        // 计算 L1 范数 (绝对误差的和)
        sumL1 += error;

        // 计算 L2 范数 (误差的平方和)
        sumL2 += error * error;

        // 计算 ∞ 范数 (最大误差)
        Linf = std::max(Linf, error);
    }

    // 计算 L1 和 L2 范数
    L1 = sumL1;
    L2 = std::sqrt(sumL2);
}

int main() {
    // 读取 JSON 文件
    std::ifstream inputFile("input/config.json");
    json config;
    inputFile >> config;

    // 解析输入参数
    std::vector<int> nValues = config["n"]; // 这里n是一个数组
    double xmin = config["xmin"];
    double xmax = config["xmax"];
    double ymin = config["ymin"];
    double ymax = config["ymax"];
    std::string bcType = config["boundary_conditions"]["type"];

    // 解析边界条件函数
    std::vector<std::function<double(double, double)>> bcConditions;
    auto bcConditionsArray = config["boundary_conditions"]["conditions"];
    for (const auto& bcExprStr : bcConditionsArray) {
        // 假设 parseExpression 函数返回一个可调用的函数对象
        auto bcCondition = parseExpression(bcExprStr);
        bcConditions.push_back(bcCondition);
    }

    // 解析源项函数
    std::string sourceExprStr = config["source_term"];
    auto f = parseExpression(sourceExprStr);

    // 解析真实函数
    std::string realFunctionStr = config["real_function"];
    auto realFunction = parseExpression(realFunctionStr);

    // 存储误差范数
    std::vector<double> L1Norms, L2Norms, LinfNorms;

    // 遍历每个 n 值，进行求解和误差计算
    for (int n : nValues) {
        std::cout << "Processing for n = " << n << std::endl;

        // 创建网格
        Grid grid(n, xmin, xmax, ymin, ymax);

        // 创建 BoundaryConditions 对象，并传入多个条件函数
        BoundaryConditions bc(bcType, bcConditions);

        // 创建 PoissonSolver 对象
        PoissonSolver solver(grid, bc, f);

        // 求解并输出结果
        solver.solve();
        Eigen::VectorXd numericalSolution = solver.getSolution();

        // 计算误差的范数
        double L1, L2, Linf;
        computeNorms(numericalSolution, realFunction, grid, L1, L2, Linf);

        // 存储误差范数
        L1Norms.push_back(L1);
        L2Norms.push_back(L2);
        LinfNorms.push_back(Linf);

        // 输出误差的范数
        std::cout << "For n = " << n << std::endl;
        std::cout << "L1-norm: " << L1 << std::endl;
        std::cout << "L2-norm: " << L2 << std::endl;
        std::cout << "∞-norm: " << Linf << std::endl;

        // 打开CSV文件进行写入
        std::ofstream outputFile("output/solution_" + std::to_string(n) + ".csv");

        // 检查文件是否成功打开
        if (outputFile.is_open()) {
            // 将Eigen向量的每个元素按行写入CSV文件
            for (int i = 0; i < numericalSolution.size(); ++i) {
                outputFile << numericalSolution[i];
                // 在元素间添加逗号，但每行最后一个元素不加逗号
                if ((i + 1) % (n+2) != 0) {
                    outputFile << ",";
                } else {
                    outputFile << std::endl;  // 换行
                }
            }

            // 关闭文件
            outputFile.close();
        } else {
            std::cerr << "无法打开文件进行写入" << std::endl;
        }
    }

    // 计算收敛率（基于 L1 范数）
    if (nValues.size() >= 2) {
        for (size_t i = 1; i < nValues.size(); ++i) {
            double n1 = nValues[i - 1];
            double n2 = nValues[i];
            double L1_1 = L1Norms[i - 1];
            double L1_2 = L1Norms[i];

            double convergenceRate = std::log(L1_2 / L1_1) / std::log(n2 / n1);
            std::cout << "Convergence rate between n = " << n1 << " and n = " << n2 << ": " << convergenceRate << std::endl;
        }
    }

    return 0;
}
